function [] = AWQobs_uncertainty(sitenm, modvarnm, obsdata, ...
              obstime, obsdepth, pos)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS FUNCTION RETURNS A DIAGRAM COMPARING SETAS RESULTS (FROM THE SETAS 
% 2-YEAR DATASET) WITH AWQ OBERVATIONS FOR A SINGLE AWQ SITE. THE DIAGRAM  
% ALSO PROVIDES AN INDIRECT UNCERTAINTY ESTIMATES OF AWQ DATA USING 
% THE STANDARD DEVIATION OF SETAS RESULTS DATA OVER PAST 6DAYS AT DATE
% OF AWQ SAMPLING.
%
% INPUT : - sitenm      : name of the AWQ site (string)
%         - modvarnm    : model variable name chosen for comparison(string)
%                         Availablable options are 'temp' or 'salt'
%         - obsdata     : matrix of AWQ weekly data (depth x date) for a
%                         single site                      
%         - obstime     : array of AWQ sampling dates 
%         - obsdepth    : matrix of AWQ sampling depth 
%         - depth       : string, either 's' for surface or 'b' for
%                         bottom AWQ data
%
% EXAMPLE: load TestmatAWQ.mat
%          AWQobs_uncertainty('B5','temp',testmatAWQ.Tempts,...
%          testmatAWQ.datesetas, 's')
% Author: Benedicte Pasquer, IMOS/eMII (http://imos.org.au/)
% email: benedicte.pasquer@utas.edu.au
% May 2013   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% DIRECTORY WHERE THE SETAS OUTPUT FILES ARE LOCATED
currentdir =pwd  ;
setasfile = dir(fullfile(currentdir,'setas_simple*.nc'));
if isempty(setasfile)
    error('No SETAS file found');
end
%INITIALIZE ARRAY OF SETAS RESULTS AND TIME
setasdata = []; setastime =[];
setas_std = [];
awqdate = []; awqdata = [];
% LOOP THROUGH SETAS MONTHLY FILES (25) TO EXTRACT THE FULL TIMESERIE AND 
% ONE STANDARD DEVIATION  VALUE PER MONTH FOR A SINGLE AWQ SITE 

for n = 1 : length(setasfile) 
    
% SOME AWQ SITE ARE LOCATED OUT OF THE BOUNDS OF SETAS GRID. NOTE THAT
% THIS IS ONLY VALID FOR SETAS_SIMPLE FILE, NOT FOR SETAS HIGH
% RESOLUTION MODEL RESULTS

  outofgrid ={'NN';'U16_17';'U19'};
  if any(strcmp(sitenm,outofgrid)) 
       
        setasdata = NaN;
        setastime = NaN;
        sitestd = NaN;
    
  else

% GET SITE LAT/LON 
        [slon,slat] = get_infoAWQ('pos',sitenm);
        if strcmp(sitenm,'CB')==1 %SITE CB MOVED AWAY FROM SHORELINE
            slon = slon + 0.004;
        end
      
% MAIN FUNCTION: GET STANDARD DEVIATION VECTOR FROM SETAS OUTPUT AT THIS LOCATION(MEAN
% OVER FIRST2 LEVELS) AND SETAS DATE VECTOR
        setasf = fullfile(currentdir,setasfile(n).name);
        [sitestd,setas_res,setas_time,dumobs,dumdate] =...
        get_setas4uncertainty(slon,slat,obstime,obsdata,obsdepth,pos,...
        modvarnm,setasf);
    
  end
  
  setasdata = [setasdata ;setas_res]; 
  setastime = [setastime ;setas_time'];
  setas_std = [setas_std; sitestd];
  awqdata   = [awqdata;dumobs];
  awqdate   = [awqdate;dumdate];
  
end

switch modvarnm
    case 'temp'
        varname = 'WaterTemp';
    case 'salt'
        varname = 'Sal';
end



% CLEAN DATA: GET RID OF NAN IN STD DATA; CONVERT DATA=0 INTO NAN 
setasdata(setasdata == 0) = NaN;
setas_std(isnan(setas_std))=[];
 awqdata (isnan(awqdata))=[];
 awqdate(isnan(awqdate))=[];
% FIGURE
figure('position',[1 1 1200 550],'Color',[.9 .98 .98])
hl = line(setastime,setasdata,'color',[.7 .7 .7]);
hd = line (awqdate,awqdata,...
    'Marker','o',...
    'Color','r',...
    'LineStyle','none',...
    'MarkerFacecolor','r');
set(gca,'Fontsize',14,'Fontname' ,'Times','Fontweigh','bold',...
    'XLim',[734382 735142],... 
    'Color', [.9 .98 .98],'YColor', [0 0 .7],'Xcolor', [0 0 .7], ...       
    'Fontname','arial','Fontsize',14,'Fontweight','bold',...
    'Position' , [0.10 0.19 0.8 0.70],...
    'Visible','on','Box','on');
hold on
% ERROR BARS
%if length(tobs_s)==length(setas_std)
    he = errorbar(awqdate,awqdata,setas_std);
    set(he,'Linestyle','none','Color','k','Linewidth',1)
%end
hold off

 


% TITLE AND LABELS
[longname4labl,unit] = get_longname(varname);
strtitre = [char(longname4labl), ' observations at site ', sitenm]; 
titre{1} = strtitre;
titre{2} = ['and CSIRO SETAS modelled surface ', char(longname4labl)];
titre{3} = [];
title(titre, ...
    'Fontname','arial','Fontsize',14,'Fontweight','bold',...
    'color',[0 0 .7]);
xlabel('Sampling date',...
    'Fontname','arial','Fontsize',14,'Fontweight','bold',...        
    'Color', [0 0 .7])
labl = [longname4labl{1} ' [' unit{1} ']']; 
ylabel(labl,...
    'Fontname','arial','Fontsize',13,'Fontweight','bold',...
    'Color', [0 0 .7])

lh = legend('CSIRO SETAS model results (hourly) ','AWQ surface observation (monthly)',...
    'Location','NorthEastOutside');

set(lh,'TextColor',[0 0 .7],...
    'Fontname','arial','Fontsize',8,'Fontweight','bold',...
    'box','off','color','none')
datetick('x','mmm-yyyy','keepticks','keeplimits')
hh = rotateXLabels(gca(),45);

set (hh,'Color',[0 0 .7])
a1=gca;
ax2 = axes('Position',[0.70 0.6 .2 0.2],...
    'XAxisLocation','bottom',...
    'YAxisLocation','left',...
    'Color','none',...
    'XLim',[0 1],'YLim',[0 1],...
    'XColor',[.9 .98 .98],'YColor',[.9 .98 .98],...
    'XTick',[ ],'Ytick',[ ] ,...
    'box','off'); 

% ADDITIONAL INFORMATION ABOUT ERRORBARS
errbrtxt{1} = 'errorbar corresponds to SETAS';
errbrtxt{2} ='standard deviation over past 6 days';
h = text(0,1,errbrtxt,'Color',[0 0 .7],'Fontname','arial','Fontsize',8,'Fontweight','bold');

set(gcf, 'InvertHardCopy', 'off');

figname = strcat(varname,'uncertainty',sitenm);
saveas(gcf,figname,'png')
close(gcf)
